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The Khinchin theorem provides the condition that a stationary pro- 
cess is ergodic, in terms of the behavior of the corresponding correla- 
tion function. Many physical systems are governed by non-stationary 
processes in which correlation functions exhibit aging. We classify 
the ergodic behavior of such systems and provide a generalization 
of Khinchin's theorem. Our work quantifies deviations from ergod- 
icity in terms of aging correlation functions. Using the framework 
of the fractional Fokker-Planck equation we obtain a simple analyt- 
ical expression for the two-time correlation function of the particle 
displacement in a general binding potential, revealing universality 
in the sense that the binding potential only enters into the prefac- 
tor through the first two moments of the corresponding Boltzmann 
distribution. We discuss applications to experimental data from sys- 
tems exhibiting anomalous dynamics. 

Anomalous diffusion | Fractional Fokker-Planck Equation | Correlation functions 
| Weak ergodicity breaking | Khinchin theorem 

The dynamics of physical observables is commonly quan- 
tified in terms of time dependent correlation functions; 
these are pair products of dynamic observables 0{t\) and 
0(t2), that are averaged over an ensemble, C(i2,ii) = 
{0{ti)0{t\)) . When such a correlation function describes a 
stationary process, C(i2,ii) is a function of the time differ- 
ence only, that is, C(t2,t\) = C(\t2 — ti\). For such pro- 
cesses the correlation function contains the dynamical infor- 
mation on other physical observables via fundamental the- 
orems [T]. For instance, the Wiener-Khinchin theorem re- 
lates the power spectrum to the correlation function; or the 
fluctuation-dissipation theorem connects correlation functions 
to linear response functions. Another well known example is 
Khinchin's theorem [2], that provides a criterion for ergodicity 
of a process in terms of the corresponding stationary corre- 
lation function. However, the stationarity property is found 
to be violated in numerous systems (see below). These sys- 
tems exhibit aging properties, that are intimately connected 
to the non-stationary behavior. Correlation functions in aging 
systems behave very differently from their stationary coun- 
terparts [3] |4] [5] [6]. For instance, such non-stationary sys- 
tems may be characterized by correlation functions of the type 
C{t2,t\) = h(ti/t2) (if t2 > t±), i.e., the two times enter as 
a quotient rather than their difference. This non-stationary 
behavior is also connected to a breaking of ergodicity in the 
sense that long time averages differ from ensemble averages of 
the same quantities [3j [7] [8] [9] [10] . The relation between cor- 
relation functions and ergodicity breaking can be quantified 
by the Edwards- Anderson parameter [3], see below. 

We here present a generalization of the Khinchin theorem 
for aging systems, relating fluctuations of time averages to the 
corresponding aging correlation function, that we calculate for 
a class of important stochastic processes. In particular we de- 
rive an analytical expression for the two-time position correla- 
tion function in the presence of an external binding potential 
U(x) based on the the general framework of the fractional 
Fokker-Planck equation [111 1121 [13] . The latter describes sys- 
tems in which the mean squared displacement in free space 
scales as (x 2 (t)} oc t a , where the anomalous diffusion exponent 
a ranges in the interval < a < 1 [121 1131 114| . Our results 



for the correlation function are quite general and could also 
be applied to continuous time random walk (CTRW) systems 
in confined geometries and mean field descriptions of models 
such as the quenched trap model [14] . In particular, we show 
that for sufficiently long times the correlation function behaves 
universally, and the dependence on the potential U(x) enters 
only in the prefactor through the first two moments of the 
corresponding Boltzmann distribution. We also discuss aging 
properties of the correlation function and the time averaged 
mean squared displacement. Moreover we demonstrate agree- 
ment of our results with experimental data. 

Physical systems displaying non-stationary behavior like 
aging and ergodicity breaking traditionally included glassy 
systems such as spin glasses [3], colloidal glasses [15], gels 
[16] . turbulent systems [17], or tracer dispersion in subsurface 
hydrology [18] , among others. More recently advanced single- 
molecule experiments reveal other types of complex systems 
with similar behavior. These are systems exhibiting anoma- 
lous diffusion and slow, non-exponential relaxation dynamics 
[141 1121 [T3l 119] . For instance, they include blinking quantum 
dots [20] [21] [22j [23], or biologically relevant systems. The lat- 
ter include subdiffusion of tracer particles in living biological 
cells [21 [251 EE EH [28] [29] or reconstituted crowding sys- 
tems 30, 31, 32, 33 , protein conformational dynamics [34] . 
or the motion of bacteria in a biofilm [35] . Here we show how 
the knowledge of the aging correlation function allows us to 
quantify the non-ergodic behavior of the process. 

Correlation functions and ergodicity. The concept of ergod- 
icity plays a central role in statistical physics. Boltzmann's 
ergodic hypothesis relates the time average of a physical ob- 
servable O to its ensemble average and states that both are 
equal in the long time limit: 



0=]im 



0(t')dt' = (O). 



[1] 



where (•) represent an ensemble average with respect to a 
steady state probability distribution. The connection between 
the ergodic hypothesis and correlation functions was estab- 
lished by Khinchin [2] , a work of considerable impact to sta- 
tistical physics [30]. Khinchin's theorem states that an ob- 
servable 0(t) is ergodic if the associated correlation function 
is "irreversible", in the sense that if 0(t) fulfills 



lim (o(h + A)o(ti)) = (oy 



(A = |t 2 — ti| 



[2] 
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then Eq. JT] holds. In the derivation of Khinchin's theo- 
rem it is assumed that the process is stationary and the sys- 
tem reached a steady state, C(t2, ii) = (0(t2)0(ti)) = C(A). 
Thus the irreversibility condition corresponds to the require- 
ment C(A) = (O) 2 for A — > oo. Note that for stationary 
systems [37] irreversibility is a broader concept than the er- 
godic hypothesis, and that therefore Khinchin's theorem can- 
not always be reversed. The case of power-law decaying corre- 
lation functions was treated previously based on the general- 
ized Langevin equation [38 , 39 . This gives rise to a stationary 
process, such that Khinchin's theorem holds for this kind of 
anomalous dynamics. 

Having in mind the aforementioned non-stationary sys- 
tems, the following question arises: Does irreversibility im- 
ply ergodicity also for aging systems? And if not, what 
theorem replaces Khinchin's? To answer these questions we 
quantify ergodicity in terms of second moments, namely if 

{^O 2 ^ = (O) 2 the dynamics is ergodic, as originally pointed 

out by Khinchin and (O) = (O) always holds for processes 
where (0(i)) is constant for long enough time. From Eq. [lj 
we find for the second moment 
1 

T 2 



o(ty 



(0(t 2 )0(ti))dt 2 dti 



3] 



Assuming aging behavior for the correlation function in the 
above form (0(t2)0(ti)) = h{t\/t2) (t2 > ii) (i.e., we use the 
aging regime as the "steady state" of the system similarly to 
invoking stationarity in Khinchin's theorem) Eq. [3[ becomes 

2 

7 1 



o{t) 



h(t 1 /t 2 )dt 1 dt 2 . [4] 

'o 

After the substitution t\ + 1| = r 2 and ti/t2 = tan(f?) we find 

[5] 



(O 2 ) = J\{z)dz. 



This is a general relation between the fluctuation of the time 
averaged and ensemble averaged correlation functions. For 
ergodicity to hold, we consequently require that in the long 
time limit the condition 



(O) 2 = / h(z)dz 
Jo 



[61 



is fulfilled. Simultaneously we can rewrite the irreversibility 
condition IimA->oo h(t/t + A) 

{Of 



(O) 2 in the form 



lim h(z) 
z->0 v ' 



7] 



Thus, to fulfill Eq. J6] the condition J7] is not sufficient, and 
Khinchin's theorem does not hold for aging systems. Namely 
irreversibility does not imply ergodicity in aging system. We 
define a generalization of Khinchin's theorem stating that in 
the case of irreversible dynamics the condition for ergodicity 
in an aging system is given by 



lim h(z) 

z— >o 



h(z) dz. 



[8] 



The knowledge of the correlation function resolves the ques- 
tion of ergodicity for a system. More importantly, by help of 
Eq. [5]] the correlation function allows one to quantify the 
fluctuations of the time averages and ergodicity breaking: 



o ) - (oy = / [h(z) - h{o)\ dz, 



[91 



where we assumed irreversibility. We now turn to the calcula- 
tion of the correlation function in an aging system and prove 
that it is irreversible. The result will be shown to exhibit uni- 
versal features which, together with Eq. ]5J, imply a generic 
behavior of the fluctuations of time averages. 



Model 

Fractional dynamics. Anomalous diffusion in an external po- 
tential U(x) — — J x F(x')dx' is governed by the fractional 
Fokker-Planck equation [11] IT2] 



df{x,t) 
dt 



= K a0 D t a L F p f(x,t), 



[101 



in which the Fokker-Planck operator Lpp 



J3_ 

dx 



SM 1 

_k B T mp \ 



+ 



-g^z includes drift and diffusion terms [40] . K a is the anoma- 
lous diffusion coefficient of dimension m 2 /sec", and fcsT mp the 
thermal energy. The fractional Riemann-Liouville operator 



Dl~ a z(t) 







z(t') 



dtr(a) J (t-t'y 



-dt', (o < a < 1) [11] 



involves long-range memory effects [121 I13| . Eq. \ \W\ \ de- 
scribes the time evolution of the single-time probability den- 
sity function f(x, i) and has been studied extensively for dif- 
ferent potentials [121 1131 HTj . Eq. [11 IJ can be derived as the 
long-time limit of a continuous time random walk model, in 
which the local probability to jump left and right is biased by 
the external potential U(x) [421143] . 

The fractional Fokker-Planck equation | ll()| | can be 
rephrased in terms of the Langevin equations [441 1451 1461 147] 



dx(s 
ds 



K 



F(x) + V (s) (a), 



dt(s) 
ds 



= <j(8) (b), [12] 



where s is an internal time (unit-less) and t is the physi- 
cal (laboratory) time. In Eq. [1121] the noise rj(s) is white 
and Gaussian with zero mean (r)(s)} = and autocorrelation 
(r](s)r)(s )) = 2KS(s — s ); K > is the diffusivity for the nor- 
mal diffusion process in internal time s. Conversely the noise 
oj(s) represents an asymmetric Levy-stable process of order a 
such that the probability density function of s is [411 1451 146] 



1 (K, 
Pt{s) = - 
a 



t 



l/a 

K j 



Ka 

K 



l/rx 



,1/a 



[131 



L a is a one-sided Levy stable probability density function with 
Laplace transform exp( — X a ) [51]. The representation of the 
fractional dynamics by means of the coupled Langevin equa- 
tions | I12I | is usually termed subordination [44]. With Eq. ([l"3]) 
it is easy to show [41] that statistical properties of x(t) are in- 
dependent of the choice of if (one may set K — 1 for example). 
Equations of the form | I12I | are useful simulation tools [461147] . 
and were used to investigate multiple-time probability density 
functions 45,48,49 and correlation functions 50, .511 . 



Two-point probability density functions. The relation between 
the solution of the fractional Fokker-Planck equation and its 
Brownian counterpart via subordination can be used to de- 
rive the m-point probability density functions for a subdiffu- 
sion process [45ll48[[49], In particular, the 2-point probability 
density function is given by [45] 

f(x2,t2;Xl,tl) — 



dsi / ds2n(s2,t2;s 1 ,t 1 )f M (x2,S2;x 1 ,s 1 ), [14] 
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where n(s2,t 2 ; si,ti) represents the 2-point probability den- 
sity function of the inverse Levy-stable process s(t), and 
fu{x2,S2',Xi,s\) is the 2-point probability density func- 
tion of the corresponding Markovian process x{s) defined in 
Eq. | |12| | (a). The Laplace transform of n is given by [45] 



n(s2, A 2 ; si, Ai) 

"{(£)«- 



exp (- Sl A^ 
A1A2 

ai) [A" - A" + AJ] 



■) e( S2 - Sl )Af [A Q - A2] e - A ? (32 - si) ^} 



exp 



+- 



Ai A2 



(Jf) 2 8(8i - s 2 )A? [A Q - A" ] e~ Af < S1 - S2 '^ 



[151 



where A = Ai + A2 and 8 is a step function. Here and in the 
following, the Laplace transform for the variable pairs ti <H> Ai 
and t2 <-> A2 is denoted by the explicit dependence on the re- 
spective variable. 

Knowing that x(s) is a Markov process, the 2-point prob- 
ability density function fM(x 2 ,S2',xi,si) of the process x(s) 
is obtained from the well known property /m(x2, S2; xi, si) = 
Pm(x2, S2\xi, si)Pm(xi, si\xo,0), where xo is the initial posi- 
tion of the particle and Pm {xi, si\xo, 0) is single time proba- 
bility function for the process x(s) found by solving Eq. H2f aV 
We decompose this expression in the form 

fM(X2,S2',XX,Sl) 

= 6(s 2 - si)Pm(x2, S2 - si\xi,0)Pm(xi,si\xo, 0) 

+8(si - s 2 )P M (xi,s 1 - s 2 \x2,0)Pm{x2,s 2 \x ,0). [16] 

From Eqs. | |14| |. |)15fl . and JTBJ we obtain the subdiffusive 
probability density function in terms of the Markovian 2-point 
probability density function Pm in Laplace space via 



f(x 2 , A 2 ; x 1 ,\ 1 )X 1 X 2 K a /K 



K 



in Laplace space. Laplace inversion (see Ref. [52] for details) 
delivers the final result for the 2-time correlation function 



{X{t2)x{tx)) ~ [(X 2 ) B - (x)i 



B(ti/t 2 ,a, 1 - q) 
T(a)r(l - a) 



+ (x)%, 



[19] 

valid for t2 > ti ^> (1/ K^fii) 1 ^ , /ii being the smallest non- 
zero eigenvalue of the Fokker-Planck operator Lfp [ID], and 



B(z,a,b) = 



V 



i (i-y) b 1 dy 



\20] 



is the incomplete Beta function [53]. The symbol (■) b denotes 
an average over the Boltzmann distribution 



{x") B = 



x n e -U{*)/k B T dx 



e -V{*)/h B T dx _ 



[211 



Fig. [T] shows the sigmoidal behavior of the position autocor- 
relation function. It is important to note that Eq. (|19|l im- 
plies that the process is irreversible since in the limit A — > 00 
Bit ^ + A, a, 1 - a) -t and (x(t 2 )x(ti)) -» (x)%. 

The result [1191] has some remarkable properties. Thus the 
external potential U (x) enters only via the prefactor, and only 
through the first two moments of the corresponding Boltz- 
mann distribution. The long time behavior of {x{t2)x{t\)) is 
universal and depends only on the ratio t 2 /t\. 

We note that Beta function behavior for a correlation func- 
tion was found previously for a simple two-state renewal model 
with power-law sojourn times on both states 3 20 22 23 , 54 
155] . While our process is clearly not a two state process, the 
universal behavior of expression | |19| | is due to the separation 
between the physical process in space and the associated tem- 
poral process. Such a separation is exactly the idea behind 
the subordination of time, Eq. [1121] (b). The temporal process 
t(s) yields the time evolution governed by the waiting times 
between successive jumps. Due to the assumption of annealed 
disorder it is independent of the current particle position. It 
converges as a function of the number of jumps (oc s) due 
to the generalized central limit theorem, corresponding to the 
long-time limit s — > 00 in Eq. | |12| |(b1. The limiting behav- 
ior of t(s) is therefore a Levy stable law which is underlying 



= [A" - A a + \%] Pm(xi, -^A^xq, 0)5(^2 - xi) 
+^{A2 [A a - A?] Pm(x 2 , ^-\%\ Xl ,0)P M ( Xl , -?-A a \x , 0) 



+A? [A Q - A" ] P M (xi, ^\?\x2,0)Pm(x 2 , J-A> ,0)}. [17] 

Here Pm (x, -^-X a \xo, 0) denotes the Laplace transform 
(^X a ^ s) of P M (x, s\x , 0), and similarly for A a . Eq. pTI ] 
is the final result connecting the 2-point probability density 
function / in Laplace space to its Markovian counterpart. 



Results 

Position autocorrelation. Eq. [ 1171] allows one to calculate 
general two-point correlation functions for subdiffusive sys- 
tems governed by the fractional Fokker-Planck equation [1101] . 
In particular, the position-position correlation function is 




(x(X 2 )x(X 1 )) 



xiX2f(x 2 ,X 2 ;xi,X 1 )dx 1 dx2 [18] 



Fig. 1. Position correlation correlation function C(t2,t\) = (x(t2)x(t\) from 
Eq. II19II (solid line) as function of the scaled time A/t for a = 0.525, where 
A = |t2 — ti [. The values of the Boltzmann moments are (x 2 ) g = 0.217 and 
(x) 2 B = 0.102. We also show data (o symbols) for the magnetization M(t+ A, t) 
from the thermoremanent magnetization experiment reported in Ref. 57 from mea- 
surement in a Cuo.94Mno.06 sample, showing aging behavior. The shown data 
correspond to the longest measured waiting time, t = 10000 sec. 
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the subdiffusion dynamics. Conversely, the spatial process 
explores the external potential and is not affected by the dis- 
order if observed as a function of the internal time s. In fact, 
as function of s the process corresponds to normal diffusion in 
an external field, and so the process converges to Boltzmann 
statistics characterised by the binding properties of the exter- 
nal potential. We note that the result | |19| | for the correlation 
function mirrors the convergence of both temporal and spatial 
processes, and is independent of the microscopic properties of 
the model (e.g. the shape of U{x)). To obtain the behavior 
when one of the processes has not converged one needs to use 
the full correlation function with a non-trivial time depen- 
dence [52] . The latter depends on all eigen- values of Lfp and 
is cumbersome. 

Properties of the two-time position correlation. The correla- 
tion function | |19| | displays a number of noteworthy features: 

(i) Aging behavior 

The correlation function C(t + A,t) = (x(t + A)x(t)) ex- 
hibits aging since its time dependence is of the form A/t. 
Aging behavior was indeed observed in many complex systems 
[3l ll6|[57] . for instance, in thermoremanent magnetization ex- 
periments 56, 57], in which the measured relaxation of the 
magnetization M(t + A, t) is proportional to the spin correla- 
tion function, according to generalized fluctuation-dissipation 
theorems [581 155[ [31. Accepting our stochastic theory as an 
approximation for the spin system behavior we used Eq. | I19I | 
to fit the thermoremanent magnetization data from Ref. [57] . 
The result of the fit is presented in Fig. [1] We observe good 
agreement between the data and our Beta function results 
over the entire time window, with a slight discrepancy at short 
times. We note that the use of a non-zero value for the fit- 
ting parameter (x) B in Eq. | |19| | for the zero external field 
relaxation of the magnetization is consistent with observed 
asymmetric magnetic fluctuations in thermoremanent magne- 
tization experiments [55] as opposed to the naively expected 
zero average behavior. Fitting with the Beta function for the 
measured correlation function does not necessarily yield in- 
sight into the physics of the system, but classification of aging 
with particular fitting functions might be profitable step (as 
the well known functions, such as Cole-Cole plots, are useful 
in the classification of dielectric relaxation). 

(ii) Time-averaged position 

We quantify ergodicity, or the departure from ergodicity, of 
the system by measuring fluctuations of the time averaged 
position 

W) = ~ [ x(t')dt'. [22] 
* Jo 

Combining Eqs. [5]] and [1191] we see that 

lim (W) 2 ) = (1 - a)«z 2 >B - (x)%) + (x)l, [23] 
where we used the relation 

f B(z,a,l-a)dz={l-a)T(a)r{l-a). [24] 
Jo 

This result was previously obtained for a CTRW process [60] . 
Clearly \x(i) y 7^ {x}% when a < 1, thus we observe weakly 
non-ergodic behavior. In the limit a = 1 ergodicity is restored. 



(iii) Edwards- Anderson parameter 

The Edwards-Anderson parameter was previously used to 
quantify the degree of weak ergodicity breaking in the con- 
text of spin-glasses [3]. It is defined as 

q E A = lim lim C{t + A,t). [25] 

A— »oo t— >oo 

In our current framework, for the case of a symmetric poten- 
tial the Edwards- Anderson parameter becomes 

qEA = { 0, a > 1 ' [26] 

reflecting irreversibility of our process which is still non- 
ergodic. Conversely, interchanging the limits we find that 

lim lim C(t + A,t) = 0. [27] 

t— >oo A— >oo 

reflecting the aging character of the system. Note that 
Eq. | |26| | indicates that qEA is determined by the Boltzmann 
distribution and is independent of a in the nonergodic phase. 

(iv) Time-averaged mean squared displacement 

From Eq. | |19| | we also obtain the behavior of the ensem- 
ble average of the time averaged mean squared displacement 
[611 162]. Namely, from a time series x(t) recorded in single 
particle tracking experiments one can define the time aver- 
aged mean squared displacement 

5HA) = ¥ ^£ A [x(t + A)-x(t)] 2 dt, [28] 

where T is the overall measurement time. At finite mea- 
surement time T even in the Brownian limit the quantity 
S 2 (A) is a random quantity depending on the particular tra- 
jectory. Performing an additional ensemble average, for a 
Brownian system the role of the lag time A in the long mea- 
surement time limit is completely interchangeable with the 




10 ' icf 2 10° 10 2 10 4 10 6 

A 



Fig. 2. Simulated behavior of the time-averaged mean squared displacement 
^(5^(A)^ for a continuous time random walk process with waiting time distribu- 
tion tp(t) ~ r a /£ 1_t " a , in an harmonic binding potential U (x) = X 2 (□), and in 
a box with reflecting boundary conditions and size 2 (o). The anomalous diffusion ex- 
ponent is a = 1/2, and the measurement time was chosen as T = 10^ (a.u.). We 
also chose k B T mp = 0.1, and K 1 / 2 = 0.0892. Without fit, the lines show the 
analytic results for the transition from the initial linear lag time dependence ~ A , 
Eq. [f32ll f— and ), to the long lag time behavior ~ A 1_a , Eq. [T3TT (—). 

In both cases (x) g = 0. At long lag times ^5 2 (A)^ /(x 2 )b exhibits universal 
behavior, independent of the external field. 
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regular f-dependence in the corresponding ensemble averaged 
mean squared displacement, for example when U(x) = 
liniT->oo S 2 (A, T) — 2KA. In presence of a confining poten- 
tial one would naively expect the mean squared displacement 
to saturate, as observed for a Brownian system. However, 
evaluating the ensemble average of <5 2 (A) 

+ A) - x(t)y^dt 



5 2 (A) 



T- A 



' [(x 2 (t + A)) + (x 2 (t)) - 2C(t + A, t)] 
T — A 



[291 



we find the a priori surprising result : For regular diffusion in 
a binding potential one obtains a saturation for long times, 
as for anomalous motion in the case of the ensemble average 
[121113] . In contrast, for the time averaged mean squared dis- 
placement in our anomalous system from Eq. J29] we find 
from the Beta function expansion 



B(t/(t + A), a, 1- 

r(a)r(i - a) 

the behavior 



1 - 



sin(a7r) 
(1 — a)ir 



(r5 2 (A)) ~ (V ) fl - (x)%) 



2 sin(a7r) 
(1 — a)a>TT 



[301 



[311 



valid in the limit (A/T) < 1 and for A > (l/if a ^i) 1/a [52]. 
Instead of the naively expected saturation, the time averaged 
mean squared displacement grows as a power with exponent 
(1 — a). Only when the lag time approaches the measurement 
time T this power-law growth stops, and the function dips to 
the ensemble averaged value. We note that the A 1 " scaling 
was recently reported for the case of a particle in a box [63] . 

(v) Numerical analysis of position autocorrelation. 

Fig. [2] shows, over a large time span, the time-averaged 
mean squared displacement 



^<5 2 (A)y of a subdiffusing particle 

in (i) an harmonic potential, and (ii) in a box with reflecting 
boundaries. The initial particle position was chosen to be at 
the bottom of the potential and in the center of the box, re- 
spectively. At short lag times A we observe the linear scaling 



<5 2 (A) 



2K a A 



r(l + 



[321 



with the lag time A. In this result only the dependence on 
the overall measurement time T bears witness to the fact 
that the underlying stochastic process is subdiffusive. Seem- 
ingly paradox, the lag time A enters linearly, in contrast 
to the associated ensemble averaged mean squared displace- 
ment (x 2 (t)) ~ 2K a t a /r(l + a). However, this is the result 
of the weak ergodicity breaking of the process, as shown in 
Refs. [611 162] . The free particle behavior at short A is an ex- 
pected result, which can be obtained from scaling arguments 
or explicitly from the full correlation function [55]: at suf- 
ficiently short times the particle does not yet feel the con- 
finement due to the reflecting boundaries, or it does not yet 
experience the restoring force of the potential, respectively. 

The ^r5 2 (A)^> ~ A regime holds for scales of the lag time A 

that fulfill K a A a <C L 2 in the example of the box, where L 
is the size of the box. For a general confining potential, the 
turnover time is non-universal and is dependent on all non- 
zero eigenvalues fj, n of the Fokker-Planck operator Lfp |40| . 
Thus, at times A ^> (1/ 'K a ^\Y^ a a transition occurs to the 




Fig. 3. Single trajectories for the motion in the harmonic binding potential 
U(x) = X 2 , for the same parameters as in Fig. [2] (thin lines). The symbols 

(□) represent the simulations result for ^<5 2 (A)^. The scatter between individ- 
ual trajectories is distinct and resembles qualitatively typical experimental results 



5 2 {A)j ~ A 1_Q regime, Eq. | |31| |. We stress again that, in 
contrast to normal diffusion, no saturation is found at long lag 
times, and <r5 2 (A) > continues to grow for any A < T. Only 



as A approaches to the measurement time T, we obtain the 

convergence ^r5 2 (A)^> ->■ (a; 2 ) + (a; 2 (0)). 

In Fig. [3] we show the simulations result for a number of 
individual trajectories in an harmonic binding potential, dis- 
playing significant scatter. This scatter between individual 
trajectories is, again, a result of the weak ergodicity breaking 
of the underlying stochastic process, and resembles qualita- 
tively the behavior observed in single particle tracking experi- 
ments [2U [25] [261 EH IM ES]. The amplitude of the scatter can 
be quantified in terms of the dimensionless random variable 

£ = S 2 (A)/ ^<5 2 (A)^, the relative scatter of the time average 

with respect to its ensemble mean. Using arguments similar 
to [62] we can show [52] that the distribution of £ is given in 
terms of a one-sided Levy stable law in the form [62] 



lim <pa(£) = TTTT, 'a 



r(i + a 



l/a 



[331 



where the stable density l a {t) is defined in terms of its Laplace 
image exp(— u a ). Note that the random variable £ is in the 
denominator of l a , and therefore the associated width is finite. 
For instance, for the case a = 1/2 used in Figs. [2]and[3]we find 
the Gaussian law Hult-s-oo <f>i/2(Q = (2/7r) exp (— £ 2 /2) . Fi- 
nally, in the Brownian limit a = 1, the distribution converges 
to the sharp behavior 0i(£) = 5(£ — 1), restoring ergodicity 
in the sense that no more scatter occurs. The distribution of 
£, given by Eq. (|33fl . is the same for both unbounded and 
bounded anomalous diffusion. This is simply due to the men- 
tioned fact that temporal and spatial process are uncoupled. 



Discussion 

Correlation functions are a standard tool to experimentally 
probe the temporal evolution of a system. They provide in- 
formation on how the present value of a physical observable 
influences its value in the future, and are therefore impor- 
tant indicators to the specific process that governs the sys- 
tem's dynamics. The significance of the correlation function 
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behavior for the fundamental concepts in physics is revealed 
through Khinchin's theorem which provides a condition for a 
stationary process to be ergodic in terms of the correspond- 
ing correlation function. Herein we derived a generalization of 
Khinchin's theorem for a class of non- stationary aging process. 
We provide not only a generalized condition for ergodicity for 
such processes in terms of the corresponding aging correla- 
tion function but also quantify the deviations from ergodic- 
ity. For the broad class of non- stationary processes described 
by Eq. [1101] we derived analytically the time dependence of 
two-point correlation functions for subdiffusing particles un- 
der situations of confinement. In particular we revealed a 
universal behavior for the two-time position correlation func- 
tion involving the incomplete Beta function. All features of 
the confining potential enter the correlation function solely 
through the prefactor in terms of the first and second mo- 
ments of the associated Boltzmann distribution. Of course, 
the expression for the correlation function in Eq. [1191] is not 
restricted to a position correlation function, but can be used 



to describe a very general class of quantities, e.g., a potential 
energy correlation function. The generality of our results is a 
direct consequence of the convergence of the Markovian pro- 
cess x(s) in the jump space, and the ubiquitous role of Levy 
statistics due to the generalized central limit theorem. 

We note that the behavior found for the time-averaged 
mean squared displacement in binding potential (or confined 
geometry) recovered herein could be verified in a single par- 
ticle tracking experiment specifically for micro-movement in- 
side a yeast cell, in which the particle was tracked indirectly 
in an optical tweezers setup, that eventually exerts a Hookean 
restoring force on the tracer particle [64]. Similarly anoma- 
lous diffusion of macro molecules in a cell is always confined 
by the cell boundary, indicating that our theory of confined 
anomalous diffusion may have wide applications. 
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